function [mg fx fy] = im_gradient(src_img)
    function dgx = gen_dgauss_x(sigma, haftsize)
        dgx_size = 2 * haftsize + 1;
        dgx = zeros(dgx_size, dgx_size);
        for m = 1:dgx_size
            for n = 1:dgx_size
                u = m - haftsize -1;
                v = n - haftsize - 1;
                dgx(m,n) = (-u/(2*pi*sigma^4)) * exp(-(u^2 + v^2)/(2*sigma^2));                
            end
        end                
    end
sigma = 1;
hsize = ceil(sigma*3);
hx = gen_dgauss_x(sigma, hsize);
hy = hx';
fx = imfilter(src_img, hx,'replicate');
fy = imfilter(src_img, hy,'replicate');
mg = sqrt(fx.^2 + fy.^2);
end